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Abstract 

A new Runge-Kutta-Nystrom method, with phase-lag of order infinity, for 
the integration of second-order periodic initial-value problems is developed 
in this paper. The new method is based on the Dormand and Prince Runge- 
Kutta-Nystrom method of algebraic order four[l|. Numerical illustrations 
indicate that the new method is much more efficient than the classical one. 
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1. Introduction 

In this paper we study a special Runge-Kutta-Nystrom method of Dor- 
mand et a/.[l| for integrating systems of ODEs of the form 

-^ = f{tMt)) (1) 

for which it is known in advantage that their solution is periodic or oscillating. 
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Several authors in their papers {for example see [3,7-10]) have developed 
Runge-Kutta-Nystrom methods with the purpose of making the phase-lag of 
the method smaller. 

The phase-lag of a method, first defined by Brusa and Nigro [2,] at 1980. 
Van der Houwen and Sommeijer [sl proposed second-order m-stage methods 
(with m = 4,5, 6) and phase-lag order g = 6, 8, 10 respectively. They also de- 
rived some third-order methods with phase-lag order 6, 8, 10. In js], Chawla 
and Rao have constructed Numerov-type methods with minimal phase-lag 
for the numerical integration of second-order initial-value problems. Simos et 
al. [sl obtain fourth-order Rungje-Kutta-Nystrom with minimal phase-lag of 
order eigth. He also derived in [9| a Runge-Kutta-Fehlberg method of order 
infinity. 

In the present paper and based on the requirements of infinite order of 
phase-lag, we will construct a phase-fitted four-stage Runge-Kutta-Nystrom 
which is based on the coefficients of the well-known Runge-Kutta-Nystrom 
Dormand et al. [l|] method of algebraic order 4. 

2. Phase lag analysis for Runge-Kutta-Nystrom methods 

The general m-stage method for the equation 



is of the form 

(0) = „ , 



fitMt)) (2) 



i 

Un = U^r^\ Un = Un-1 + bjfj, (3) 



where 



fi = /(tn-1 + Cih, Un-1 + hCiUn-1 + X] ^) (4) 

and Ci = and = 1 

The above expressions are presented using the well-known Butcher table, 
given below: 
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Cm ^m,l (-^in,2 ■ ■ ■ ^m,m—l 

bi 62 • • • bm-1 bm 

bi 62 • • • bm 

Table 1: m-stage Runge-Kutta-Nystom method 



In order to develop the new method, we use the test equation, 



{ivfuit) =^ u"{t) = -v^u{t), V eR 



(5) 



By applying the general method ([2]) to the test equation ^ we obtain 
the numerical solution 



Un 

hUr, 



huQ J 



D 



z = vh, 



(6) 



where A,B,A',B' are polynomials in z^, completely determined by the pa- 
rameters of method ([3]) 

The exact solution of is given by 



where 



u{tn) = cri[exp(zf)]" + cr2[exp(-if)]". 



0-1,2 = ± ^^^] or ai^2 = W\exp{±ix)- 



(7) 



2^ - V 
Substituting in (JTj), we have 

u(tn) = 2\a\cos{x + ^-2) 



(8) 



Furthermore we assume that the eigenvalues of D are Qi, ^2, and the conse- 
quent eigenvectors are [1,^2]"^, 

where Vi = A' / {pi — B'),i = 1, 2. The numerical solution of (I5l) is 



ClPi + C2P2 , 



(9) 
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where 

V2U0 — huQ ViUq — huQ 

Ci = , C2 = . 

f 1 -V2 Vi- V2 

If pi,p2 are complex conjugate, then ci,2 = \c\exp{±.iw) and pi^2 = 
\p\exp{±ip). By substituting in Qj, we have 

Un = 2\c\\p\'^'cos{w + np). (10) 

From equations ([8]) and ffTOj) we take the following definition. 



Definition 1. (Phase-lag). Apply the RKN method ^ to the general 
method Then we define the phase-lag = z ~ p. If ^{z) = O^z"^^), 
then the RKN method is said to have phase-lag order q. 

In addition, the quantity a{z) = 1 — |p| is called amplification error. 
Let us denote 

R{z'') = tr{D) = A{z^) + B'{z^) 

Q{z^) = det{D) = A{z^)B\z^) - A'{z^)B{z^) (11) 
where z = vh. From Definition [1] it follows that 

^{z) = z- arcoss f — ^i£L) , |p| = y/Q{^. (12) 



2x/g( 



z^' 



We can also put forward an alternative definition for the case of infinite 
order of phase lag. 

Definition 2. (Phase-lag of order infinity). To obtain phase-lag of order 
infinity the relation = z — arccosi ^/^ )= must hold. 

3. Derivation of the new Runge-Kutta-Nystrom method 

In this section we construct a 4-stage explicit Runge-Kutta-Nystrom 
method (presented in Table 1), based on R^z"^) and Q^z"^). Now let us rewrite 
R and Q in the following form 

R{z^) = 2 - riz"^ + r2Z^ - r^z^ + ... + Viz"^' = 

Q{z^) = 1 - qiz^ + q2Z^ - qsz^ + ... + qiZ^' = (13) 
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By computing the polynomials A, B, A' , B' and therefore it! and Q in 
terms of RKN parameters we obtain the following expressions 

Alz"^) = l+64a4,3a3,2a2,i2;* + (-&4a4,2a2,i-&3a3,2a2,i-^4a4,3a3,i-&4a4,3a3,2);s^ + 
(6202,1 + 6404,1 + 6404,2 + 6303,1 + 6404,3 + b3a3^2)z'^ + (-64 - 61 - 63 - 62)2;^ 

B{z^) = 1- ha4^sa3^2C2z'^ + (6404,3C3 + 6404,202 + 6303,202)2;'^ + (-63C3 - 64 C4 - 

6202)-^^ 

A^z"^) = 6404,303,202,1^^ + (-6303,202,1 - 6404,303,1 - 6404,303,2 - 6404,202,1)^^ + 
(6202,1 + 6303,1 + 6303,2 + 6404,1 + 6404,2 + 6404,3)-^^ + (-64 - 62 - 61 - 63)2;^ 

B'{z^) = 1-6404,303,2022;^ + (6404,303 + 6404,202 + 6303,202)2;"' + (-63O3 - 64O4 - 

6202)2;^ 

R{z'^) = 2+6404,303,202,12;^ + (-6303,202,1 - 6404,303,2 - 6404,202,1 - 6404,303,202 - 
6404,303,1)2;^ + (6202,1 + 6303,2 + 6404,3 + 6303,202 + 6404,303 + 6404,202 + 6303,1 + 
6404,1 + 6404,2)2;^ + (-63 - 62 - 6303 - 6404 - 6202 - 64 - 61)2;^ 

Q{Z^) = l + (-64a4,3a3,i62C2 - 6404,202,16303 - 6202,16404,303-6303,202,16404 + 

^303,16404,202 — 6303,16404,202 — 6404,303,202,1 + 6404,202,16303-616404,303,202 + 

6404,16303,202 - 6404,16303,202 + 6404,303,202,1 + 616404,303,202+6303,202,16404 + 

64O4,3O3,l6202+62O2,l64O4,303)2;^ + (-64O4,3O3,l-63O3,2O2,l -6404,202,1-6404,303,2- 
616404,202-616303,202-636404,202+626404,303-6202,16303-6202,16404-6404,16303- 
6404,16202-6404,26303 + 6303,26404+6404,16303+6404,16202+6404,26303-626404,303 + 
6202,16404 + 6303,16404 + 6303,16202-6303,16404-6303,16202-6404,36202-6303,26404- 
646303,202-616404,303+646303,202+616404,303+6404,36202+616404,202+616303,202 + 
636404,202 + 6202,16303 + 6303,202,1 + 6404,303,1+6404,202,1+6404,303,2 — 6404,303,202) 
2;^ + ( - 64 63 O3 + 6462 O2 - 62 64 C4 - 63 64 O4 +^62 64 O4 + 63 62 O2 + 6364O4 - 61 63 O3 +64 63 O3 + 
6163C3 — 6462O2 — 6162O2 — 6164O4 + 6162O2 + 6263O3 + 6164O4 — 6362O2-6263O3 — 
6202,1 - 6303,1 - 6303,2 - 6404,1 - 6404,2 - 6404,3 + 6202,1 + 6404,1 + 6404,2 + 
6303,1+6404,3 + 6303,2 + 6303,202 + 6404,303 + 6404,202)2;^ + (-62 - 64O4 + 62 - 
64 + 61+63 - 62O2 - 61 - 63 - 63O3 + 64)2;^ 

where z = uh 
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As it has already been defined, in order to have phase-lag of order infinity, 
the following relation must hold: 



<!>{z) = z- arccos ( — ^ \^ ) = (14) 



By applying R{z^) and Q{z'^) to the formula of the direct calculation of 
the phase lag (fT2l) and substituting the following coefficients that have been 
used by Dormand et al. in [l| : 



^21 
k 



1 7 119 

32' 
1 

4' 
1 

14' 
1 

14' 



1000' 




"32 = 


500' - 14' 


"42 = , 

27 


7 

^^=10' 


C4 = 


= 1, 






27 


h 


25 
189' 


&4 = 0, 




h 32 

02 = , 

81 


h 


250 
567' 


^^ = 54' 





After satisfying relation ( fT4l) . we have: 



^(z) = z — arcoss 



5 1 

a4 3 = z ^ -(54621 -4793320/ + 99172960/ 

5292 289^4 _ 6800^2 + 40000;^^ v 

+ 5179680 (sin {z)f - 768268800 / + 4043520 z^ (sin {z)f 

+ 1866240000 - 559872000 (sin {z)f + 24 (-654383577600 / 

+ 212348252160000 - 1366377865200 / - 1710031785 

+ 89285428680 - 202307339750400 z^ (sin (2))^ 

+ 2023399802880000 / (sin (^))^ - 2015539200000000 / 

+ 581660870400 / (sin {z)f + 1319799592800 / (sin {z)f 

+ 1710031785 (sin {z)f - 89285428680 (sin {z)f 

+ 46578272400 (sin (z) )^ + 72722707200 z^ (sin (z) )^ 

- 10040912409600 / (sin (2))^ + 544195584000000 (sin {z) f 

- 7860602880000 / (sin {z) f + 6046617600000000 

- 6590813184000000 (sin (2))^)^/^) (15) 
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The Taylor expansion series for 04,3, which is given from the above formula 
is : 

25 _ 43 2 _ 1531 4_ 3273029 g 

"^'^ ~ 189 ~ 2400 ^ ~ 30240000 ^ ~ 36288000000 ^ 

59772887431 o , , 

9699782400000000 ^ ^ 

4. Numerical examples 

In this section we will apply our method to three problems. We are go- 
ing to compare our results with those derived by using the high order method 
of embedded Runge-Kutta-Nystrom 4(3)4 method of Dormand and Prince 
(see l^ij). 

One way to measure the efficiency of the method is to compute the ac- 
curacy in the decimal digits, that is —logio{maximum error through the in- 
tegration intervals) 

acc{T) = —logio{max\u{tn)—Un\), where tn = l+nh, n = l,2, 
and u{t) is the vector of the solution. 

Table 2 shows the accuracy for the two methods. In our computations we 
have two step values, for Problems 1 and 2, h = 0.025 and h = 0.050, and 
for Problems 3 and 4, h = 0.25 and h = 0.50. 

Problem l.(Inhomogeneous equation) 

= t + Z/2 _ 1 sin t , M = 1, ^'0=^^ + 1, 

dt^ 

where t > and v = 10. 

The analytical solution is u(t) = cos{ut) + sin{ut) + sin(t) 
Problem 2. (Two-Body problem) 



U = Z 



(u2 + z2)3/2' (^2 + ^2)3/2 

where u{0) = 1, u'{0) = 0, z{0) = 0, z'{0) = I and v 
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Our method Dormand and Prince method 





T=100 


T=1000 


T=5000 


T=100 


T=1000 


T=5000 


Problem 1 














h=0.025 


4.2 


3.2 


2.5 


2.3 


1.3 


0.6 


h=0.050 


2.7 


1.7 


1.0 


1.1 


0.2 


-0.3 


Problem 2 














h=0.025 


7.3 


5.9 


4.6 


6.5 


5.1 


3.8 


h=0.050 


6.0 


4.4 


3.1 


5.2 


3.6 


2.3 


Problem 3 














h=0.25 


5.7 


5.4 


5.4 


4.2 


4.1 


4.1 


h=0.50 


4.2 


3.9 


3.9 


2.9 


2.8 


2.8 


Problem 4 














h=0.25 


5.2 


4.3 


3.4 


3.5 


2.5 


1.6 


h=0.50 


3.8 


2.8 


1.9 


2.3 


1.8 


0.4 



Table 2; Accuracy for the maximum absolute error for problems 1-4 



The analytical solution is u{t) — cos{t) and z{t) = sin{t) 
Problem 3. (Duffing equation) 

^ -u{t) - {u{t)f + Bcos{ut) 
where B = 0.002 and ly = 1.01. 

The analytical solution is u{t) — Aicos{i't) + AsCos{3iyt) + A5Cos{5iyt) + 
Arcosllut) + Agcos{9i't) 

where Ai = 0.200179477536, A3 = 0.000246946143, A5 = 0.000000304014, 
Ar = 0.000000000374 and Ag = 0.000000000000 

Problem 4 -(Franco and Palacios problem) 

u{t) e C u{Q) = 1, u'{Q) = (1 - ]-e)i, 

- cos{t) + ^etsin{t) + i[sin{t) — ^etcos{t)] 



d u[t) . . ... 

_— ^ = -u{t) + eexp[it), 

where e = 0.001 and u = 1 
The analytical solution is u{t) 
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5. Conclusion 

A new fourth order Runge-Kutta-Nystrom method with phase-lag of 
order infinity is developed in the present paper. The new method is based 
on the very well known classical Dormand and Prince fourth algebraic or- 
der Runge-Kutta-Nystom method. The numerical results show that the new 
method is much more efficient for integrating second-order equations with 
periodic oscillating behavior than the classical one. 
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